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Analytic methods show stabihty of the stationary accretion of test fluids but they are inconclusive in the 
case of self-gravitating stationary flows. We investigate numerically stability of those stationary flows onto 
compact objects that are transsonic and rich in gas. In all studied examples solutions appear stable. Numerical 
investigation suggests also that the analogy between sonic and event horizons holds for small perturbations of 
compact support but fails in the case of finite perturbations. 
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I. INTRODUCTION 

Investigation of the spherical accretion onto compact ob- 
jects starts with seminal works of Hoyle, Lyttelton, and Bondi 
m [3[] describing the infall of dust matter onto a surface 
of a star moving through the interstellar medium. The first 
hydrodynamical analysis of spherical accretion was presented 
in 1952 by Bondi |4] who considered a spherically symmet- 
ric flow of a polytropic perfect fluid in a Keplerian gravita- 
tional potential. A general-relativistic model of accretion in 
the Schwarzschild space-time has been developed by Michel 
ill]. In all these works the test fluid approximation has been 
adopted. 

The first general-relativistic model with self-gravitating 
steady fluids has been analyzed by Malec [6]. This work 
has been later continued by Karkowski, Kinasiewicz, Mach, 
Malec, and Swierczyiiski |7], and resulted in finding a whole 
family of steady transsonic solutions of Einstein equations de- 
scribing a self-gravitating cloud of gas accreting onto a central 
compact object (a black hole in this particular case). The most 
striking fact about these solutions is that, given fixed asymp- 
totic parameters of the model (such as the size of the cloud, 
speed of sound at the outer boundary, and the total asymp- 
totic mass of the system), there exist two different transsonic 
solutions corresponding to the same accretion rate — one, for 
which most of the mass is contained in the central object, and 
the other, where the amount of mass contained in the accreting 
fluid constitute almost all the mass of the entire configuration. 
The first class of the solutions contains a subset of the test 
fluid flows found by Michel. Other solutions are new. Let 
us stress, that although these solutions are not available in a 
closed form, many of their parameters and properties can be 
inferred by analytical means. Such flows might be associated 
with the Thorne-Zy tkow stars [8,, ^ or quasistars ifToll . 

We should also mention here the work of Papadopoulos 
and Font [11], who constructed a general-relativistic hydro- 
dynamical code capable of simulating self-gravitating flows 
and applied it to the strongly perturbed Michel's solution. The 
fluid perturbation was treated as self-gravitating but the back- 
ground solution corresponded to the test fluid regime. A sim- 
ilar investigation, but in the context of radiation hydrodynam- 
ics, has been done by Zampieri et al. fl^]. They investigated 
the stability of solutions in Schwarzschild space-time found 
by Nobilietal. |13I] . 

The first proof of the stability of the transsonic accretion in 



the Newtonian and relativistic, spherically symmetric cases, 
given by Moncrief |14], was restricted to the test fluid ap- 
proximation. Analytic methods are inconclusive in the case 
of self-gravitating accretion flows ifTsll . Below we report the 
results of a numerical analysis of the stability of both branches 
of solutions found in |7]. An interesting by-product of this in- 
vestigation is that in the nonlinear regime sonic horizons are 
movable and the signal can get out from within the sonic hori- 
zon in the original steady flow. That hints to the limited valid- 
ity of the formal analogy between sonic and event horizons. 

A stability analysis of Newtonian accretion solutions with 
self-gravitating flows has been done for axially symmetric 
perturbations. In all examined cases we have observed the 
stable behavior of both aforementioned branches of solutions: 
those corresponding to the test fluid as well as those where the 
mass of the fluid dominates over the mass of a central object. 

The order of the forthcoming sections is as follows. Section 
II gives a short description of the two branches of accreting 
flows in spherically symmetric space-times. Section III intro- 
duces the dynamical equations of motion in the form adopted 
in the numerical code. Section VI reports results concern- 
ing stability of accreting flows. The next section shows that 
sonic horizons can be penetrated from within by large per- 
turbations. This suggests the limited validity of the analogy 
between sonic and event horizons. The stability of Newto- 
nian accretion under axisymmetric perturbations is shown in 
section VIII. Obtained results are briefly reviewed in section 
VIII. 



II. STEADY SOLUTIONS 

A general spherically symmetric space-time can be de- 
scribed by the line element 



-N^df ■ 



■ adr^ + [dO^ + sin^ Od^'^) 



(1) 



where A^, a, and R are functions of the coordinate radius r and 
the asymptotic time variable t. 

The extrinsic curvature of a slice of constant time t of the 
space-time with the metric given by ([TJ has the following non- 
zero elements: K'^ = dra/{2aN), ^ ^ djRURN). Ac- 
cordingly the trace of the extrinsic curvature can be written as 
\xK = N-^^J\n[^[^R^). 

Similar calculations can be performed for the two-spheres 
of constant radius r embedded in a given temporal slice. The 



2 



1.6 10"" 
1.4- 10-" 
1.2- 10-''' 
1.0- 10-" 
■S 8.0- 10-'* 
6.0- 10-'* 
4.0- 10-'* 
2.0- 10-'* 
0.0- 10° 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

'"lluid/'"iot 



FIG. 1: The dependence of the accretion rate m on the mnnid. Here 
mtot = 1- 



result for the trace of the extrinsic curvature (twice the mean 
curvature of the surface) is A: = 2d,-R/ (r Vq')- 

In what follows we will consider the evolution of a spher- 
ical cloud of perfect fluid accreting onto a central object and 
described by the energy-momentum tensor 



r'"' = (p+p)MV + pg' 



(2) 



where p is the pressure, p the energy density, and the four- 
velocity of the fluid. 

We start working with the comoving gauge, so that - 
— u!^ - Q (there exist a suitable geometric condition im- 
posed on extrinsic curvatures Kj for such a choice of coordi- 
nates [@]), and introduce a function U = RKg = djRjN. It 
has the meaning of the spatial part of the fluid four-velocity 
computed in the reference frame (?', r') that has been obtained 
by the transformation (f, r) i-> {f - t,r' - R{t, r)). We also 
introduce the quasilocal mass, which can be easily expressed 
as 



m(R) 



•47r 



f R' 

Jr 



pdR'. 



Here mtot is the total asymptotic mass of the configuration, and 
Rco denotes the size of the accretion cloud. Another important 
quantity, the local speed of sound a is defined by 

Symbols x and k are used here to denote derivatives 



which have to be computed according to the assumed equation 
of state; the quantity h = {p+p)/n is the specific enthalpy. For 
a barotropic equation of state the above definition reduces to 
= dpi dp. 

We will search for the solutions of the Einstein equations 
according to the following notion of stationarity. The accre- 
tion rate m - dfin, computed at a given areal radius R should 



be constant in time. Similarly, other hydrodynamical quanti- 
ties like the fluid velocity U, pressure p, energy density p, etc. 
should satisfy d,'X = 0, where X - U,p,p,... These assump- 
tions can be satisfied only approximately as the accreting fluid 
contributes to the growing mass of the central object and the 
whole configuration must change in time. We will, however, 
see that they lead to solutions characterized by very small ac- 
cretion rates both in the test and in the heavy fluid regime. 
Thus, for a very long time (much larger than a characteris- 
tic dynamical time understood as a time required by a sound 
wave to travel across the cloud) the motion of the accreting 
gas remains almost unchanged, and the above assumptions are 
justified a posteriori. 

Under these assumptions the set of Einstein partial differ- 
ential equations and the equations expressing the conservation 
of the energy-momentum tensor reduces to a set of ordinary 
ones, namely 



d , IN 
— In 
dR 



\kRl 



\6n A 



Bn 

p + p' 



kR^2 



2m 

1 + U^. 

R 



(3) 



Here A and B are for integration constants. The baryonic 
density n is defined as a function assuring that Vij(nu^) - 0. 
When deriving these equations we have tacitly assumed that 
the equation of state is of a barotropic form p = p(p) or, equiv- 
alently, p = pin). In the following we will specialize to the 
polytropic equation of state p - Ktf . 

In addition to these equations a suitable set of boundary 
conditions has to be specified. As such we usually choose the 
size of the accretion cloud R^o, its total mass mtot, the asymp- 
totic value of the baryonic density «oo (alternatively poo can be 
also used), and the asymptotic value of the speed of sound Ooo. 

Such boundary conditions still provide a family of the so- 
lutions parametrized by the value of J/m, or the constant A. 
In this work we are only interested in the so-called transsonic 
solutions — those for which far away from the central object 
the fluid is subsonic while in the central parts of the cloud it 
falls supersonically. Such a solution passes through the so- 
called sonic point, defined as a location where 2U/(kR) = a. 

The analysis presented in |Ql is concerned with accretion 
flows onto a central black hole which will naturally appear in 
the model if we continue to integrate the equations from the 
outer boundary inwards until the apparent horizon is passed 
by. In the terms of this paper the apparent horizon is defined as 
a surface on which the optical scalar 6+ - kR/2 + U vanishes. 
According to this definition, the radius of the apparent horizon 
Rbh and the mass iubh - '«(^bh) (the mass of the black hole) 
satisfy the standard relation Rbh - 2mBH- Apart from the 
mass of the black hole we also define the fluid mass as mguid = 
mtot - niBH- 

Fig.[T] shows the dependence of the accretion rate on the 
ratio of mouid/mtot for the specific case of a sequence of poly- 
tropic models with Roo = 10*, mtot = 1 and Ooo = 0.1. (The ra- 
tio of mfluid/mtot scales linearily with «„.) Clearly, the slowly 
accreting regime corresponds either to a situation with a very 
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small mass of the fluid as compared to the mass of the cen- 
tral black hole or to converse case where almost entire mass 
is in the form of fluid. This can actually be proved analyti- 
cally [T*]. Moreover, it can be shown that the maximum of in 

corresponds to OTfluid/'«tot =1/3- 

The attempts to investigate the stability of the new branch 
of massive solutions analytically have failed to give any con- 
clusive results (see, e.g., 1,15.1 ). Due to this fact we devote this 
paper to study this numericaUy. 



III. DESCRIPTION OF THE NUMERICAL CODE 

The dynamical code used to investigate the stability of 
steady solutions was constructed in a similar fashion to the 
one described in 1 16]. It is a modern version of a High Resolu- 
tion Shock Capturing (HRSC) scheme based on the Godunov- 
type methods developed for solving the equations of hydrody- 
namics. 

For the construction of the code the polar gauge has been 
used. The metric is assumed to be of the form 

ds^ = -a'^dt^ + X^dR^ + (d0^ + sin^ Od^"^) , (4) 

where the lapse a and X are functions of the areal radius R 
and time t. Following ITtIi we introduce the Lorentz factor 
W = au' and the three- velocity v' - u'/W. The equations of 
conservation of the energy-momentum tensor V^T'"' - and 
the baryonic density Vftinu^) = can be now written as 



The numerical scheme is stabilized by a suitable choice of 



d,q + j^dR [aXR'-F] = - (d, \nX) q. 
Here q denotes a vector of conserved quantities 

q = (D,S,Tf = (nW, hIiW^vr, nWihW - 1) - pf , 
F stands for the flux vector 

F = (nWv*, nhW^VRV^ + p, nW{hW - l)v*)^ , 
and E denotes the source terms 







Time derivatives of conserved quantities q are computed 
according to the following version of the method of lines 



(dq 

I dt 



iaXR^F)M/2 - {aXR^F)i-i/2 



XiR^ARi 



+ iaZ-—q 



(5) 



where lower indices refer to spatial cells (shells of constant 
radius). Values of q corresponding to the subsequent time- 
step are obtained using standard Runge-Kutta methods. 



the numerical fluxes F 



/+I/2- 



In most of the modern HRSC 



schemes numerical fluxes at the cells interfaces are computed 
based on the solutions to the local Riemann problems that 
arise naturally between each of the cells' interfaces. In the 
following qL,i+i/2 and qR,,+i/2 will denote left and right Rie- 
mann states at the i+ 1/2 interface. In order to provide higher 
order of the spatial accuracy the states qL,/+i/2 and qR,/+i/2 are 
computed as follows 

= q/ + S;(/?/+ 1/2 
q!+l/2 - Q'+l + S'+l - Ri+i) ■ 

Here Ri and /?,+i/2 are the positions of the cells' centers and 
interfaces respectively. The slope limiters S , are defined by 



S, = minmod 



q/+i - q/ q, - q, i 



[Rni - Ri Ri - Ri-i I 
where the "minmod" function has been introduced as in ifTsIl 

{a if \a\ < \b\, ab > 0, 
b if \a\ > \b\, ab > 0, 
iffl/7<0. 

As a method to compute numerical fluxes F,+i/2 we have 
used two versions of a scheme proposed originally by Donat 
and Maraquina in Uftl . both of them being based on the spec- 
tral decomposition of the Jacobian dF jdq consisting of three 
eigenvalues Ap, left eigenvectors 1^ and right ones r,,. In the 
original version of the algorithm one starts by computing the 
following variables 



and 



/ - l'' r. / - l'^ r. 



=l(;-F(qL), < = l'^F(qR), 



where p numbers the eigenvectors of dF jdq. Now, for each 
p, the signs of A'^ and A'^ are inspected. If both eigenvalues A'^ 
and A'^ are positive, we define 

while for both A'^ and A'^ having negative values, we set 

If the signs of the two eigenvalues A'^ and A'^ differ, one defines 

l^''lmax(L,P) = max(|/l(^|, and 

f_ = i« + |^nmax(L,pH). 

The total numerical flux is now computed according to the 
formula 

F,..,/2 = 2K< + ^'^r^R)- 
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Apart from this algorithm, known in the hterature as the 
Maraquina's flux formula, we have also implemented its slight 
modification described in \.20i . Here the numerical fluxes are 
computed as 



F,+i/2 



i^F(qL)+F(qR) + 



We have not observed any significant difference in the code 
performance between these two schemes. 

A spectral decomposition of the Jacobian dF/dq can be eas- 
ily obtained analytically. Its eigenvalues read 



Xv'^ + a 



Xv'^-a 



X + a\'R X - avR 

Let us introduce the following quantities 

1 - VflV« 



•,2' 



1 - Vr/1± ' 



and 



With their help the right eigenvectors of 5F / 5q can be written 

as 

r =l^v 1-— f 

r± = (l,hWX^A^:H^,hWJl:t - if ■ 
For the left eigenvectors we have 
W 



•7C- 



-[h-W,Wv'^,-W) , 



' hWJl:^ (vfi - X^A^) -VR+ 'KX^J[^A^ ' 
1 - 'KJi^, 
-vr + 'Km^X^A^ 



In the spherically symmetric case and the polar gauge the 
Einstein equations can be reduced to just two ordinary dif- 
ferential equations in R, that can be subsequently solved by 
quadratures provided that the hydrodynamical equations are 
known (see e.g. ll2lll '). The first of these equations 

dRin = AnR^ (nhW^ ~ p) " 4nR^{D + t) 

gives the radial derivative of the quasilocal mass, related to 
the metric function X by 



X = 



1 



(6) 



2m 



The second one provides the radial derivative of the logarithm 
of the lapse 

dR\na = + ^ttTJ {nhW^VRV^ + p)^ 

These equations are integrated (numerically) at each time- 
step, that is after new values of the hydrodynamic quantities 
have been obtained. Notice that the equation for mass m has to 
be solved before we attempt to integrate the equation for In a. 
For a closed system one can fix the total mass of the system 
and integrate the first equation starting from the outer bound- 
ary Roa- Then the second equation can be integrated in the 
same way assuming, for instance, that the lapse at the outer 
boundary is given by the standard expression known from the 
Schwarzschild solution in the polar gauge. 



a(R^)^ /l- 



2m{Roo) 



Ra 



Let us remark, however, that such assumption for the lapse 
is not necessary for the proper behavior of the hydrodynamic 
part of the code. 

The Einstein equations yield also the following expression 

d,X = -AnnhaW^XvRR, 

which is required in order to establish the source terms needed 
by the evolution scheme Q. 

Recovery of the primitive quantities (n, v*, p) from the con- 
served ones {D, S,t) is performed every time-step by means 
of the Newton-Raphson technique (we solve an equation for 
the pressure p). 



IV. CODE TESTS 

Our numerical code has been tested on the spherical shock 
reflection problem ll22ll . the Michel solution for spherical ac- 
cretion in the Schwarzschild space-time [Sj], and models of 
spherical poly tropic stars ll23ll . 

The first test checks the validity of the hydrodynamical part 
of the code in spherical symmetry. Initial data for this problem 
consist of a spherically symmetric flow of perfect fluid with 
a constant, negative radial velocity vq. The initial baryonic 
density distribution is also constant and the internal energy 
density is set to a negligibly small value. Such initial data 
evolve by producing a strong shock wave appearing at r = 
and propagating outward with velocity equal to 

(r-i)Woivoi 



Wo+l 

where Wo is the Lorentz factor corresponding to vq. 

The second test checks the validity of the implementation 
in a case of the test fluid accretion occurring in the fixed 
Schwarzschild background. The initial data for this test can 
be easily obtained by solving an algebraic equation for each 
value of the areal radius. 
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The initial data for the third test, namely a static solution 
describing a polytropic star, have to be computed by solving 
the Tolman-Oppenheimer-Volkov equations for a polytropic 
equation of state |23]. In this case satisfactory results have 
been obtained usin g a R unge-Kutta scheme of 8-th order by 
Hairer and Wanner ll24ll . 

All these tests have been passed as desired, convincing us 
that all parts of the code work properly. 



V. INITIAL CONDITIONS 

In order to construct the initial data for the main study of 
this paper we have taken numerical, transsonic solutions to the 
equations (O and added a perturbation in velocity. The way 
of obtaining such solutions is simple but not entirely obvious; 
we will describe it briefly here. 

Instead of solving the equations in R we introduce a new 
independent variable ( - I /R, which results in a grid that 
becomes naturally dense in the inner regions of the cloud and 
relatively coarse outside. Next we introduce the following set 
of dependent variables 



r'(t, r) - R(t, r). This leads to the line element 



1 = dR'R'^p, 
Jr 

Jr k^R' 



yi - 167r 

^3 = a^-, 



and express the equations (O in terms of y\, y2, ^3, and 
This yields a set of three differential-algebraic equations that 
can be integrated starting from - ^IRoo towards increasing 
values of f (that is from the outer boundary to the center of 
the cloud) provided that the values of and A are 

specified. 

We solved these equations using Daspk — a solver for a sys- 
tem of differential and algebraic equations developed by Pet- 
zold. Brown, Hindermarsh, and Li [25i, 126ll . 

Usually, the solution found in this way will not pass through 
the sonic point. We can, however, search for the transsonic so- 
lution by exploiting the fact that it can be integrated to values 
of ^ corresponding to the region inside the apparent horizon. 
(This is not true for other solutions which break down out- 
side the horizon.) Thus, in order to find a transsonic flow, 
we have implemented a bisection method, which looks for a 
value of A giving the maximal ( at which the corresponding 
solution crashes. After the proper value of A had been ob- 
tained, we could confirm that the appropriate solution indeed 
passes through the sonic point. 

Such a solution is expressed in the coordinate system where 
the time coordinate is the comoving time. Before treating it as 
a possible initial data one has to express the fluid velocity in 
the polar gauge. By adopting the areal radius as the radial 
coordinate we have changed from the comoving gauge metric 
given by the line element ([T]i to the coordinates f'(f, r) - t. 



2^ / 2 \2 

dt'^ -2UN\ — \ dt'dR + 
kR 



+ j dR^ + R^ [deP- + sin^ ed(p^) . 



Here the four-velocity of the fluid reads 



U, m" 



0. 



The transition to the polar coordinates with the metric of the 
form (HI can also be easily done. The function X is given by 
^ and the velocity - u'^/W — which is one of the dynami- 
cal variables in our code — can be written as 



U 



In this way we have obtained four functions m(R), v*(^), n(R), 
and e{R) that correspond to stationary flow. 



VI. STABILITY OF ACCRETING SELF-GRAVITATING 
FLOWS 

The code described in the preceding section has been used 
in order to evolve perturbed steady accretion flows. 

In addition to the initial data one has to specify suitable 
boundary conditions. At the inner boundary outflow condi- 
tions have been assumed (the fluid was allowed to fall in- 
ward). This boundary was positioned outside the apparent 
horizon but always in the supersonic zone of the accretion 
cloud, so that the outflow conditions could be easily imple- 
mented. Thus, our numerical models resemble also a situation 
in which the accreting fluid is glued to the solid surface of the 
central body. Such an approach has been adopted for instance 
in 1I27I1 . where the processes of the radiation transport through 
the accreting medium are also taken into account. 

The outer boundary was kept fixed using the values ob- 
tained from the initial solution (the ghost zones were filled 
with the appropriate initial values). This last condition cannot 
be easily relaxed. After setting the outer boundary condition 
to an outflowing one, even the test fluid solution can be unsta- 
ble. 

As the first step in the stability analysis we performed a 
consistency check, by assuming initial data for the dynamical 
equations to be equal to m{R), v'^(R), n{R), and e{R) inherited 
from the steady flow solution. It appeared that the evolution 
did not produce any noticeable changes; this confirms the va- 
lidity of the assumption of stationarity. Next, we introduce 
an additional perturbation. In our case it was applied to the 
velocity (the other three initial data m(R), n{R), and e(R) 
are the same as in the steady flow); \^ was perturbed by a 
bell-shape profile of a sine wave restricted to one half of its 
period. Such perturbation, initially located outside the sonic 
radius, produces two signals: one traveling outwards and one 
towards the center. If the initial amplitude of the perturbation 
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FIG. 2: Evolution of small perturbations. Here iiq denotes the density 
in the unperturbed flow. The snapshots show the density contrast 
(w - no) /no in the chronological order. 



is very small and its support is relatively narrow, the signal 
propagating inwards passes through the sonic point and even- 
tually disappears through the inner boundary (it falls onto the 
central body). 

Fig.|2]shows the evolution of tiny perturbations of compact 
support. Here, in order to visualize any changes, we were 
forced to plot the contrast of density, i.e., (« - «())/«()> where 
«o refers to the unperturbed flow. In our case this value is of 
order 10"^, demonstrating the quality of the method by being 
able to handle solutions up to such precision. This, however, 
requires some fine tuning in the post-processing of the data. 
It is, for instance, well known that the very process of inter- 
polating the initial solution onto a new grid introduces small 
numerical errors, which can easily be observed at high preci- 
sion. Thus, in order to eliminate such effects, we have evolved 
the unperturbed flow to get the values of no used to calculate 
the density contrast at different times. 

The background solution used in this example coiTe- 
sponded to the following parameters: The exponent in the 
polytropic equation of state was set to F = 1 .4, the asymptotic 
mass of the whole configuration has been normalized to unity, 
the outer boundary of the cloud was placed at Roc = 10'', the 
asymptotic baryonic density was equal «oo = 1.6 ■ lO"'*^, and 
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FIG. 3: Evolution of the perturbed density. The snapshots are placed 
in the chronological order. The profile corresponding to the unper- 
turbed flow is depicted with a dotted line. 



the asymptotic sound velocity was set to =0.1 (here and in 
all other results we have adopted the gravitational units with 
c - G = 1). For such a solution the sonic point and the appar- 
ent horizon were located at Rt = 0.82 and R^h = 0.34 respec- 
tively, and the mass contained in the fluid constituted the bulk 
of the entire mass, namely mfluid = 0.83. The accretion rate 
for this solution reaches a very small value of m = 2.6 ■ 10"'**, 
thus the growth of the central object can be neglected during 
entire simulation (this fact has been confirmed independently 
by allowing the central mass to grow according to the actual 
value of ;«). 



VII. SONIC HORIZONS VERSUS APPARENT HORIZONS 

For large initial perturbations the situation can be more 
complex, as will be discussed below. A discontinuous solution 
with shocks can develop, and we can observe some reflection 
of the signal that was initially propagating inwards. 

Fig.|3] shows snapshots from the evolution of the pertur- 
bations applied to the same background solution as before. 
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FIG. 4: Reflection of the incoming signal from the inner parts of the 
accretion cloud. The velocity Xv" is drawn with the solid line, while 
the dotted line depicts values of -a. Intersections of the two graphs 
correspond to the sonic points (see the discussion in text). 



In each of the graphs on Fig. [3] the perturbed density profile 
has been plotted over the profile corresponding to the steady 
solution. One can clearly see the stable behavior on these 
plots even though they are limited to the small radius range 
of R < 0.5 ■ 10^ and relatively short evolution times. The 
original simulations have been performed during much longer 
times confirming the stability of the flow. At late stages of the 
evolution the initial perturbation was indistinguishable from 
the numerical noise. 

All these numerical results suggest the stability of the ac- 
cretion flow also in the regime where the mass of the fluid is 
large. Naturally, numerical simulations of this kind can never 
replace a strict mathematical proof, because we cannot inves- 
tigate the whole family of possible initial data (meaning both 
steady solutions and perturbation profiles). 

The reflection of the strong incoming signal from the in- 
ner parts of the accreting cloud is somewhat surprising. For 
a small and compact perturbation entering the supersonic re- 
gion in the center of the cloud it should not be possible to be 
reflected and reach the subsonic region again — this would re- 
quire the perturbation to travel with a velocity which, relative 



to the unperturbed flow, should be greater than the speed of 
sound. Due to that fact, in analogy with the black hole hori- 
zon, the term "sonic horizon" has been coined to name the 
surface bounding the supersonic region. 

A careful inspection reveals that, contrary to the event hori- 
zons suiTounding black holes, this notion can only be approx- 
imate. For strong nonlinear perturbations the concept of an 
unperturbed background solution loses meaning. The values 
of the local sound speed a change due to the perturbation and 
new sonic points can appear. In addition, the speed of a strong 
shock wave is no longer limited to the local speed of sound. 

Such behavior has been illustrated on Fig.H] Here we are 
still dealing with the same initial solution, perturbed slightly 
more to show the whole phenomenon in a clearer way. The 
quantity Xv* is plotted with the solid line, while for the graph 
of -a a dotted line has been used. A part of the initial per- 
turbation develops into the shock propagating inwards which, 
after some time of evolution, creates additional sonic points 
(i.e., intersections of graphs of Xv* and -a) lying outside the 
original "sonic horizon." When the perturbation reaches the 
inner parts of the cloud it even destroys the "original" sonic 
point, which is being replaced by the newly created one. The 
reflected signal really becomes visible in the broad structure 
outside the new sonic point, or to say this differently, the out- 
ermost sonic point moves toward the center releasing the out- 
going perturbation, which then propagates freely outwards. 

The behavior of small perturbations of compact support 
confirms the standard interpretation of sonic horizons as 
analogs of event horizons in the linear regime. However, re- 
sults concerning strong perturbations suggest that one should 
be careful in formulating the analogy (quite common in the 
so-called analogous gravity models) between the sonic hori- 
zon and the event horizon [28, 29]. 

In summary, our numerical results reported in the last 
two sections suggest that in both cases — of small and large 
perturbations — steady accretion is stable. The amplitudes of 
the perturbations decrease, and after sufficient time we are left 
with the background solution. This can be demonstrated for 
both accretion regimes: solutions with a large mass in the cen- 
ter and having little fluid and systems with large amount of gas 
but possessing light compact cores. 



VIII. STABILITY OF THE NEWTONIAN ACCRETION 

In this Section we shall report results obtained with the use 
of a version of the Prometheus code by Fryxell, Miiller, and 
Arnett ifs^ . adapted to our purposes. Prometheus is a Newto- 
nian HRSC hydrodynamical code implementing the original 
PPM reconstruction scheme developed by Colella and Wood- 
ward isill . It has been extensively used for simulating such 
astrophysical phenomena as supernova explosions and it is ca- 
pable of simulating self-gravitating flows both in spherically 
and axially symmetric cases. As the initial data we have used 
solutions to the Newtonian equations for the steady flow that 
include self-gravitation of the accreting fluid. An equation of 
state was polytropic, p = Kp^. A point mass nip was assumed 
to exist atR-Q. 
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FIG. 5: Evolution of axially symmetric density perturbations of the 
steady accretion cloud. The plots show the spatial distribution of the 
density at times 5 • 10', 5 • lO'*, and 2.5 • 10' respectively. 



Results concerning spherically symmetric perturbations 
have been already presented in |32]. They are qualitatively 
similar to the relativistic results reported in the preceding sec- 
tions. The accreting flow was stable both in the test and fluid- 
rich regimes. Large incoming perturbations were being re- 
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FIG. 6: Continuation of the previous figure. The subsequent snap- 
shots correspond to evolution times of 4.5 • 10^, 6.5 • 10^, and 8.5 • 10' 
respectively. 



fleeted from the inner parts of the accretion cloud. We have 
also observed effects analogous to those described in the pre- 
vious chapter concerning the creation and destruction of the 
sonic horizons. Moreover, for sufficiently small perturbations 
of compact support no reflection has been noticed. 

Below we shall report studies of axially symmetric pertur- 
bations. In Prometheus simulations of two or three dimen- 
sional flows are being performed using the so-called dimen- 
sional splitting. The local Riemann problems appearing in the 
Godunov-type method are being solved separately in each of 
the dimensions, but, in order to preserve consistency of the 
method, every time-step the transversal components of veloc- 
ity are also being evolved, using the effective advection equa- 
tions. The subtle point about this procedure is that the order 
of the subsequent one dimensional sweeps is important to pro- 
vide the desired accuracy of the whole scheme. This issue has 
been discussed in detail by Strang i n ll33ll . 

The gravitational potential has to be found at each time- 
step by solving the gravitational Poisson equation. The 
Prometheus code implements a method, based on the expan- 
sion of the gravitational potential in terms of the spherical 
harmonics, that has been developed by Miiller and Steinmetz 

in. 

In our case the Euler equations of hydrodynamics and the 
Poisson equation for the gravitational field were solved on a 
spherical grid consisting of 600 zones in the radial and 180 
zones in the angular direction. 

The results presented on Fig. |5] have been obtained for the 
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following parameters of the flow. The outer boundary of the 
cloud was assumed to be located at R^o = 2-10^, the central 
mass was set to mp - 3 ■ lO-' (in the Newtonian limit gravi- 
tational units that we use here correspond to setting G = 1). 
As before, we have chosen a solution where most of the mass 
is contained in the fluid, i.e., mp/wtot = 3%. The asymp- 
totic parameters of the unperturbed flow were as follows: 
Poo = 3 ■ 10"'^, Uoo - -4 ■ 10"^, and for the parameters of the 
poly tropic equation of state we have taken F = lA, K - S-IQ'^. 

Initial data coincide with the relevant characteristics of the 
stationary flow, with the exception of the velocity. The veloc- 
ity perturbations are of the form 

6U(R 0)^1^ sinCee) for < R < Rh, 

\ 6 for R<R, or Ri, < R, 

where A is the initial amplitude of the perturbations while R„ 
and Rh describe their support. This is obviously a completely 
arbitrary choice. The only important feature of these pertur- 
bations is their nonsphericality. 

The boundary conditions have been chosen in the same way 
as in the relativistic case. At the inner boundary the outflow 
conditions were assumed, while at the outer boundary we have 
implemented "external" inflow conditions based on the back- 
ground steady flow. Fig.|5] and |6] show snapshots from the 
evolution of the density, which has been color-coded using a 
logarithmic scale. 

As expected, after a sufficient time all perturbations either 
disappear through the inner boundary or disperse outwards. 
This again suggests that the accretion flow is stable, even in 
the regime where the mass of the fluid dominates over that of 
the central object. In particular, we do not observe any frag- 
mentation of the cloud, which might be expected in the case 
of non-spherical, self-gravitating accretion flows. A simple 
theoretical argument, basing on the unproven f33\ but unrea- 
sonably effective Jeans criterion, goes as follows. The Jeans 
length can be estimated as R] = Ooo/ y/g^ ~ 6 ■ 10^. This gives 
a value three times larger than the size of the entire accreting 
cloud, i.e., Rca = 2 ■ 10*. That hints at the absence of a frag- 
mentation due to self-gravity, which agrees with our numeri- 
cal results. On the other hand it is not clear whether the Jeans 
criterion can be applied to general nonspherical perturbations 
and the instability in the general case cannot be excluded. 



IX. CONCLUSIONS 

This paper is dedicated to the discussion of the stability of 
steady accretion of perfect fluids onto compact objects, with 
an emphasis on the self-gravitation of the accreting gas. The 
recently discovered steady accretion flows are rich in the fluid 
and the backreaction effects are important. Known analytic re- 
sults do not apply to such systems. The stability investigation 
requires a proper handling of a set of nonlinear partial differ- 
ential equations and it is accessible only by means of numer- 
ical computations. Such a numerical analysis has been per- 
formed in preceding Sections using modern, high resolution 
and shock-capturing schemes, both in the general-relativistic 



case as well as in Newtonian hydrodynamics. In all examined 
cases the flows have been stable, even for large nonlinear per- 
turbations. Simulations of the evolution of large, sometimes 
even discontinuous, perturbations have led to a remarkable 
observation on the so-called "sonic horizons." In the linear 
regime, with very small perturbations, the sonic horizon can 
be viewed as a surface bounding a region from which no per- 
turbation can escape; it is impenetrable from inside, analogous 
to the event horizons in general relativity. We show that for 
strong perturbations this is no longer true. Perturbations can 
change positions of the sonic points and easily escape from a 
region initially bounded by the sonic horizon. Thus, the anal- 
ogy between the "sonic horizon" and the event horizon of a 
black hole is rather limited. 

Our analysis of the stability is restricted to the accretion 
cloud only. In all numerical simulations we had to ensure that 
the gas is being delivered to the system with a small, constant 
rate. To relax this assumption one would have to take into 
account a physical process that could be responsible for the 
feeding of the accretion cloud. One of the possibilities is to 
consider the so-called "quasistars" fl^. They consist of the 
spherical accretion cloud surrounded by a large, highly mas- 
sive stellar envelope, being a reservoir of gas necessary to sup- 
port the accretion. The other class of objects might constitute 
Thorne-Zy tkow stars JH |^ . 

In section IV we described a number of classical tests for 
numerical codes. Stationary accretion flows posess all es- 
sential elements — velocity field, selfgravitation and pressure. 
These solutions are stable. We think that they should be in- 
cluded into the standard test suite for general-relativistic hy- 
drodynamical codes. 

Spherically symmetric models of accretion must be under- 
stood as a highly idealized version of physical reality. In most 
astrophysical scenarios the inflowing gas is observed in the 
form of accretion discs deviating strongly from spherical sym- 
metry. Nevertheless, models of spherical accretion occupy an 
important place in the theoretical astrophysics as elements of 
a more complex description of astrophysical phenomena (see, 
e.g., the aforementioned works on "quasistars" by Begelman, 
Rossi, and Armitage |10]). They allow for the inspection of 
the effects caused by self-gravity in general-relativistic hydro- 
dynamics in a simple, but nontrivial, case. 

More realistic models should take into account the radia- 
tion originated and transported through the accretion cloud. 
Recently a Newtonian analysis of such processes has been 
published by Karkowski, Malec, and Roszkowski [21]. These 
works are being continued in the general-relativistic context, 
the results revealing the importance of self-gravity and its con- 
nection to the observational characteristics of the model such 
as the total luminosity of the accretion cloud or the redshift of 
the emitted radiation. 
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